;;..................................................................
;; zonal aerosol water content 
;;..................................................................
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"    

begin

  fontscale = 0.8 

;;..................................................................
;; input data 
;;..................................................................

  fna = "X0006A_0003-0007_ANN_pres.nc" 
  fnb = "T63L31K01_2010_ANN_pres.nc" 

  varn = "v300" 
  varn = "v500" 
  varn = "v700" 
  varn = "v200" 
  varn = "v850" 

  load "load.ncl" 

  fla = addfile(fna,"r") 
  flb = addfile(fnb,"r") 

  va = fla->$varn$
  vb = flb->$varn$

;;..................................................................
;; colormap 
;;..................................................................
  colormap = "WhiteBlueGreenYellowRed" 

;;..................................................................
;; open wks 
;;..................................................................
  wks = gsn_open_wks("pdf",varn)            ; open a ps file
  gsn_define_colormap(wks,colormap)      ; choose colormap

  np = 3 

;;..................................................................
;; plot configure
;;..................................................................
  plot = new(np,graphic)                       ; create a plot array
  
;;..................................................................
;; plot defination 
;;..................................................................
  res                     = True
  res@gsnDraw             = False           ; don't draw
  res@gsnFrame            = False           ; don't advance frame
  res@cnInfoLabelOn       = False           ; turn off cn info label
  res@cnFillOn            = True            ; turn on color
  res@cnLinesOn            = False              ; turn off contour lines
  res@gsnSpreadColors     = True            ; spread out color table
    
  res@gsnSpreadColorStart = 4
  ;;res@gsnSpreadColorEnd   = -3 
  
  res@lbOrientation    =   "vertical"          ; vertical label bar
  res@lbLabelFontHeightF = 0.021 * fontscale 
  res@lbTitleOn        = False
  res@lbTitleString    = "" 
  res@lbTitleFontHeightF= .015                 ; make title smaller
  ;;res@lbLabelAngleF    =    45
  ;;res@pmLabelBarParallelPosF   = 0.00 
  ;;res@pmLabelBarOrthogonalPosF = 0.12           ; move whole thing down

  res@vpWidthF             = 0.8        ; width of contour plots

  res@tmYRMode             = "Automatic"          ; turn off special labels on right axis

  res@cnLevelSelectionMode = "ExplicitLevels"   ; set explicit contour levels
  res@cnLevels            = (/0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5,0.7,1.0,1.2,1.5/) 

;;  res@cnLevelSelectionMode =  "ManualLevels"   
;;  res@cnMinLevelValF       = -10.
;;  res@cnMaxLevelValF       =  45.
;;  res@cnLevelSpacingF      =   5. 

;;..................................................................
;; Font 
;;..................................................................
  res@gsnStringFontHeightF = 0.025 * fontscale    
  res@tmXBLabelFontHeightF = 0.022 * fontscale  
  res@tmYLLabelFontHeightF = 0.022 * fontscale  
  res@tiMainFontHeightF    = 0.025 * fontscale    
  res@tiXAxisFontHeightF   = 0.025 * fontscale  
  res@tiYAxisFontHeightF   = 0.025 * fontscale  

;;..................................................................
;; individual plot 
;;..................................................................

  res@gsnRightString = "kg kg~S~-1~N~" 
  res@gsnRightString = "m s~S~-1~N~" 

  ip = 0 

  ks = 4 ;;6

  res@gsnLeftString = "Charasteristic updraft velocity" 
  res@gsnLeftString = "CAM5" 
  res@tiXAxisString = "" 
  res@tiYAxisString = "" ;;Pressure (hPa)" 
  plot(ip) = gsn_csm_contour_map_ce(wks,va(0,:,:),res) 
  ip = ip + 1 

  res@gsnLeftString = "ECHAM5-HAM2" 
  res@tiXAxisString = ""
  res@tiYAxisString = "" ;;Pressure (hPa)"
  plot(ip) = gsn_csm_contour_map_ce(wks,vb(0,:,:),res)
  ip = ip + 1


;;;  res@tiYAxisString = "" 
;;;  res@gsnLeftString = "b) HAM2 default (~F33~k~F~-Koehler 95%)" 
;;;  plot(ip) = gsn_csm_pres_hgt(wks,vzb(ks:,:,0),res) 
;;;  ip = ip + 1 
;;;
;;;  res@tiXAxisString = "Latitude"  
;;;  res@tiYAxisString = "Pressure (hPa)" 
;;;  res@gsnLeftString = "c) HAM2 with ZSR 100%" 
;;;  plot(ip) = gsn_csm_pres_hgt(wks,vzc(ks:,:,0),res) 
;;;  ip = ip + 1 
;;;
;;;  res@tiYAxisString = "" 
;;;  res@gsnLeftString = "d) HAM2 with ZSR 95%" 
;;;  plot(ip) = gsn_csm_pres_hgt(wks,vzd(ks:,:,0),res) 
;;;  ip = ip + 1 
;;;
;;;;;..................................................................
;;;;; global mean aod, wataod, burden 
;;;;;
;;;;;   cases        ttaod     wataod     burden 
;;;;; 
;;;;;   HAM1         0.140      0.105       75.2
;;;;;   HAM2         0.135      0.094       50.9
;;;;;   ZSR 100%     0.163      0.123       99.9
;;;;;   ZSR  95%     0.138      0.097       56.0 
;;;;;..................................................................
;;;
;;;  ttaod = (/0.140,0.135,0.163,0.138/)
;;;  awaod = (/0.105,0.094,0.123,0.097/)
;;;  awpct = (/74.8,69.8,75.1,69.7/)
;;;  burden= (/75.2,50.9,99.9,56.0/)
;;;
;;;;;..................................................................
;;;;; add text : global mean aod, wataod, burden 
;;;;;..................................................................
;;;  txres               = True
;;;  txres@txFontHeightF = 0.022 * fontscale 
;;;  txres@Position      = "Left"
;;;  txres@txJust        = "CenterLeft"
;;;
;;;  texta = new(np,graphic)
;;;  textb = new(np,graphic)
;;;  textc = new(np,graphic)
;;;  
;;;  do ip = 0,np-1
;;;
;;;     ;;ca = "Total AOD        " 
;;;     cb = "Water AOD      " 
;;;     cc = "Water Burden   " 
;;;
;;;     ;;printa = ca + sprintf("%3.3f",ttaod(ip)) 
;;;     printb = cb + sprintf("%3.3f",awaod(ip)) + " ("+sprintf("%3.1f",awpct(ip))+"%)"  
;;;     printc = cc + sprintf("%3.1f",burden(ip))+" Tg" 
;;;
;;;     ;;texta(ip) = gsn_add_text(wks,plot(ip),printa,-60.,250.,txres)
;;;     textb(ip) = gsn_add_text(wks,plot(ip),printb,-60.,340.,txres)
;;;     textc(ip) = gsn_add_text(wks,plot(ip),printc,-60.,290.,txres)
;;;
;;;  end do

;;..................................................................
;; create panel
;;..................................................................
  resP                  = True                   ; modify the panel plot
  resP@gsnPanelLabelBar = False 
  ;;resP@gsnPanelBottom   = 0.05
  resP@txString   = ""
  resP@gsnPanelYWhiteSpacePercent = 5
  resP@gsnPanelXWhiteSpacePercent = 0
  gsn_panel(wks,plot,(/3,1/),resP)               ; now draw as one plot


end 









